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ABSTRACT 

This paper describes a general purpose computer program for analyzing steady state and 
transient flow in a complex network. The program is capable of modeling phase changes, 
compressibility, mixture thermodynamics and external body forces such as gravity and 
centrifugal. The program’s preprocessor allows the user to interactively develop a fluid 
network simulation consisting of nodes and branches. Mass, energy and specie 
conservation equations are solved at the nodes; the momentum conservation equations are 
solved in the branches. 

The program contains subroutines for computing “real fluid” thermodynamic and 
thermophysical properties for 33 fluids. The fluids are: helium, methane, neon, nitrogen, 
carbon monoxide, oxygen, argon, carbon dioxide, fluorine, hydrogen, parahydrogen, 
water, kerosene (RP-1), isobutane, butane, deuterium, ethane, ethylene, hydrogen sulfide, 
krypton, propane, xenon, R-ll, R-12, R-22, R-32, R-123, R-124, R-125, R-134A, R- 
152A, nitrogen trifluoride and ammonia. The program also provides the options of using 
any incompressible fluid with constant density and viscosity or ideal gas. 

Seventeen different resistance/source options are provided for modeling momentum 
sources or sinks in the branches. These options include: pipe flow, flow through a 
restriction, non-circular duct, pipe flow with entrance and/or exit losses, thin sharp 
orifice, thick orifice, square edge reduction, square edge expansion, rotating annular duct, 
rotating radial duct, labyrinth seal, parallel plates, common fittings and valves, pump 
characteristics, pump power, valve with a given loss coefficient, and a Joule-Thompson 
device. 

The system of equations describing the fluid network is solved by a hybrid numerical 
method that is a combination of the Newton-Raphson and successive substitution 
methods. This paper also illustrates the application and verification of the code by 
comparison with Hardy Cross method for steady state flow and analytical solution for 
unsteady flow. 



INTRODUCTION 


A fluid flow network consists of a group of flow branches, such as pipes and ducts, that 
are joined together at a number of nodes. They can range from simple systems consisting 
of a few nodes and branches to very complex networks containing many flow branches 
simulating valves, orifices, bends, pumps and turbines. In the analysis of existing or 
proposed networks, some node pressures and temperatures are specified or known. The 
problem is to determine all unknown nodal pressures, temperatures and branch flow rates. 

This program requires that the flow network be resolved into nodes and branches. The 
program’s preprocessor allows the user to interactively develop a fluid network 
simulation consisting of nodes and branches. In each branch, the momentum equation is 
solved to obtain the flow rate in that branch. At each node, the conservation of mass, 
energy and species equations are solved to obtain the pressures at that node. 

The oldest method for systematically solving a problem consisting of steady flow in a 
pipe network is the Hardy Cross method [1]. Not only is this method suited for solutions 
generated by hand, but it has also been widely employed for use in computer generated 
solutions. But as computers allowed much larger networks to be analyzed, it became 
apparent that the convergence of the Hardy Cross method might be very slow or even fail 
to provide a solution in some cases. The main reason for this numerical difficulty is that 
the Hardy Cross method does not solve the system of equations simultaneously. It 
considers a portion of the flow network to determine the continuity and momentum 
errors. The head loss and the flow rates are corrected and then it proceeds to an adjacent 
portion of the circuit. This process is continued until the whole circuit is completed. 
This sequence of operations is repeated until the continuity and momentum errors are 
minimized. It is evident that the Hardy Cross method belongs in the category of 
successive substitution methods and it is likely that it may encounter convergence 
difficulties for large circuits. In later years, the Newton-Raphson method has been 
utilized [2] to solve large networks, and with improvements in algorithms based on the 
Newton-Raphson method, computer storage requirements are not much larger than those 
needed by the Hardy Cross method. 

The objective of the present effort is to develop: a) a robust and efficient numerical 
algorithm to solve a system of equations describing a flow network containing phase 
changes, mixing and rotation, and b) to implement the algorithm in a structured, easy-to- 
use computer program. 



MATHEMATICAL FORMULATION 


GFSSP assumes a Newtonian, steady, non-reacting and one dimensional flow in the flow 
circuit. The flow can be either laminar or turbulent, incompressible or compressible, with 
or without heat transfer, phase change and mixing. The analysis of the flow and pressure 
distribution in a complex fluid flow network requires resolution of the system into nodes 
and branches. Nodes can be either boundary nodes or internal nodes. Pressures, 
temperatures, and concentrations of fluid species are specified at the boundary nodes. At 
each internal node, scalar properties such as pressures, temperatures, enthalpies, and 
mixture concentrations are computed. The flow rates (vector properties) are computed at 
the branches. 


GOVERNING EQUATIONS 

Figure 1 displays a schematic showing adjacent nodes, their connecting branches, and the 
indexing system used by GFSSP. In order to solve for the unknown variables, mass, 
energy and fluid specie conservation equations are written for each internal node and flow 
rate equations are written for each branch. 
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Figure 1 - Schematic of GFSSP Nodes, Branches and Indexing Practice 


Mass Conservation Equation 
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Equation 1 requires that, for the transient formulation, the net mass flow from a given 
node must equate to rate of change of mass in the control volume. In the steady state 
formulation, the left hand side of the equation is zero. This implies that the total mass 
flow rate into a node is equal to the total mass flow rate out of the node. 


Momentum Conservation Equation 

The flow rate in a branch is calculated from the momentum conservation equation 
(Equation 2) which represents the balance of fluid forces acting on a given branch. A 
typical branch configuration is shown in Figure 2. Inertia, pressure, gravity, friction and 
centrifugal forces are considered in the conservation equation. In addition to these five 
forces, a source term S has been provided in the equation to input pump characteristics or 
to input power to a pump in a given branch. If a pump is located in a given branch, all 
other forces except pressure are set to zero. The source term, S, is set to zero in all 
branches without a pump. 
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Figure 2 - Schematic of a Branch Showing Gravity and Rotation 
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The two terms in the left hand side of the momentum equation represent the inertia of the 
fluid. The first one is the time dependent term and must be considered for unsteady 
calculations. The second term is significant when there is a large change in area or 
density from branch to branch. The first term in the right hand side of the momentum 
equation represents the pressure gradient in the branch. The pressures are located at the 
upstream and downstream face of a branch. The second term represents the effect of 
gravity. The gravity vector makes an angle (0) with the flow direction vector. The third 
term represents the frictional effect. Friction was modeled as a product of and the 
square of the flow rate and area. K f is a function of the fluid density in the branch and the 
nature of flow passage being modeled by the branch. The calculation of Kf for different 
types of flow passages has been described in detail later within this report. The fourth 
term in the momentum equation represents the effect of the centrifugal force. This term 
will be present only when the branch is rotating as shown in Figure 2. K rot is the factor 
representing the fluid rotation. K,. ot is unity when the fluid and the surrounding solid 
surface rotates with the same speed. This term also requires a knowledge of the distances 
between the upstream and downstream faces of the branch from the axis of rotation. 

Energy Conservation Equation 


The energy conservation equation for node i, shown in Figure 1, can be expressed 
mathematically as shown in Equation 3. 
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Equation 3 shows that for transient flow, the rate of increase of internal energy in the 
control volume is equal to the rate of energy transport into the control volume minus the 
rate of energy transport from the control volume plus the rate of work done on the fluid 
by the pressure force plus the rate of work done on the fluid by the viscous force plus the 
rate of heat transfer into the control volume. 



For a steady state situation, the energy conservation equation. Equation 3, states that the 
net energy flow from a given node must equate to zero. In other words, the total energy 
leaving a node is equal to the total energy coming into the node from neighboring nodes 
and from any external heat sources (Q,) coming into the node and work done on the fluid 
by pressure and viscous forces. The MAX operator used in Equation 3 is known as an 
upwind differencing scheme which has been extensively employed in the numerical 
solution of Navier-Stokes equations in convective heat transfer and fluid flow applications 
[3]. When the flow direction is not known, this operator allows the transport of energy 
only from its upstream neighbor. In other words, the upstream neighbor influences its 
downstream neighbor but not vice versa. The second term in the right hand side 
represents the work done on the fluid by the pressure and viscous force. The difference 
between the steady and unsteady formulation lies in the left hand side of the equation. 
For a steady state situation, the left hand side of Equation 3 is zero, where as in unsteady 
cases the left hand side of the equation must be evaluated. 


Fluid Specie Conservation Equation 

The flow network shown in Figure 1 has a fluid mixture flowing in most of the branches. 
In order to calculate the density of the mixture, the concentration of the individual fluid 
species within the branch must be determined. The concentration for the k th specie can be 
written as 
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For a transient flow, Equation 4, states that the rate of increase of the concentration of 
the k th specie in the control volume equals the rate of transport of the k th specie into the 
control volume minus the rate of transport of the k th specie out of the control volume. 

Like Equation 3, for steady state conditions, Equation 4 requires that the net mass flow of 
the k ,h specie from a given node must equate to zero. In other words, the total mass flow 
rate of the given specie into a node is equal to the total mass flow rate of the same specie 
out of that node. For steady state, the left hand side of Equation 4 is zero. For the 
unsteady formulation, the resident mass in the control volume is changing and therefore, 
the left hand side must be computed. 


Thermodynamic and Thermophysical Properties 

The momentum conservation equation, Equation 2, requires knowledge of the density and 
the viscosity of the fluid within the branch. These properties are functions of the 
temperatures, pressures and concentrations of fluid species for a mixture. Three 



thermodynamic property routines have been integrated into the program to provide the 
required fluid property data. GASP [4] provides the thermodynamic and transport 
properties for ten fluids. These fluids include Hydrogen, Oxygen, Helium, Nitrogen, 
Methane, Carbon Dioxide, Carbon Monoxide, Argon, Neon and Fluorine. WASP [5] 
provides the thermodynamic and transport properties for water and steam. GASPAK [6] 
provides thermodynamic properties for helium, methane, neon, nitrogen, carbon 
monoxide, oxygen, argon, carbon dioxide, hydrogen, parahydrogen, water, isobutane, 
butane, deuterium, ethane, ethylene, hydrogen sulfide, krypton, propane, xenon, R- 1 1 , R- 
12, R-22, R-32, R-123, R-124, R-125, R-134A, R-152A, nitrogen trifluoride and 
ammonia. For RP-1 fuel, a look up table of properties has been generated by a modified 
version of GASP. An interpolation routine has been developed to extract the required 
properties from the tabulated data. 


Friction Calculations 


It was mentioned earlier in this paper that the friction term in the momentum equation is 
expressed as a product of K f , the square of the flow rate and the flow area. Empirical 
information is necessary to estimate K f . Several options for flow passage resistance are 
listed in Table 1 . 


Table 1 - Resistance Options in GFSSP 


Option 

Type of Resistance 

Input Parameters 

Option 

Type of Resistance 

Input Parameters 

1 

Pipe Flow 

L (in), D (in), 
e/ D 

10 

Rotating Radial 
Duct 

L (in), D (in), 
N (rpm) 

2 

Flow Through 
Restriction 


11 

Labyrinth Seal 

rj (in), c (in), m 
(in), n, a 


Non-circular Duct 

a (in), b (in) 

12 

Flow Between 
Parallel Plates 

r s (in), c (in), 
L (in) 

4 

Pipe with Entrance 
and Exit Loss 

L (in), D (in), 
s/D, Kj, K e 

13 

Common Fittings 
and Valves (Two K 
Method) 

D (in), K„ K 2 

5 

Thin, Sharp Orifice 

D, (in), D 2 (in) 

14 

Pump 

Characteristics 1 


6 

Thick orifice 

L (in), D, (in), 
D 2 (in) 

15 

Pump Power 


7 

Square Reduction 

D, (in), D 2 (in) 

16 

Valve with Given 

c v 

C v , A 

8 

Square Expansion 

D, (in), D 2 (in) 

17 

Joule-Thompson 

Device 

L 0 hm’ k^, A 

9 

Rotating Annular 
Duct 

L (in), r 0 (in), 
r j (in), N (rpm) 
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1 Pump characteristics are expressed as A p = Ao + Bo"* 

A p - Pressure rise, lbf/ft 2 


m - Flow rate, lbm/sec 


COMPUTER PROGRAM 

GFSSP was developed on an IBM compatible PC using the Lahey EM32 FORTRAN 
compiler. The same source code also runs on a Silicon Graphics workstation. The code 
was developed with a modular structure to achieve two objectives. First, the code’s 
solver module was separated from the preprocessor module such that users are not 
required to write any code to develop their models. Secondly, the code can easily be 
extended to enhance its capability by adding new modules. The main routine controls all 
of the program operations and makes the decisions whether to continue or to stop the 
calculations. 

The computer program has three major parts. The first part consists of the subroutines for 
the preprocessor. The preprocessor allows the user to interactively create the flow 
network model consisting of nodes and branches. All of the input specifications, 
including the boundary conditions, are specified through the preprocessor. The second 
major part of the program consists of the subroutines that provide the initial conditions 
and then develop and solve all of the conservation equations in the flow network. The 
third part of the program consists of the thermodynamic property programs, GASP, 
WASP, and GASPAK that provide the necessary thermodynamic and thermophysical 
property data required to solve the resulting system of equations. 

Figure 3 shows GFSSP’s process flow diagram. The user runs the interactive 
preprocessor to generate the input data file. The input data file contains all the 
information necessary for the model. The solver module reads the input data file and 
produces the solution in conjunction with the thermodynamic property programs. It 
should be noted that the user interfaces with the program only through the input data file. 
The user is never required to modify the source program to develop the model. The user 
is also not required to understand the structure of the solver module in order to develop 
their model. However, for the completeness of the documentation, a flow chart showing 
the major activities of the code is described in Figure 4. 
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Figure 3. - GFSSP Process Flow Diagram 
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Figure 4. GFSSP Flowchart of Major Subroutines 



RESULTS & DISCUSSION 


The paper also describes comparison of steady state flow network simulation with the 
predictions of Hardy Cross method and unsteady flow simulation of tank blow down with 
an analytical solution. A typical water distribution network is shown in Figure 5 and 
Table 2. 
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Figure 5 - Water Distribution Network Schematic 


Table 2 - Water Distribution Network Branch Data 


Branch 

Length(inches) 

Diameter(inches) 

Roughness 

Factor 

12 

120 

6 

0.0018 

25 

2400 

6 

0.0018 

27 

2400 

5 

0.0018 

57 

1440 

4 

0.0018 

53 

120 

5 

0.0018 

56 

2400 

4 

0.0018 

64 

120 

4 

0.0018 

68 

1440 

4 

0.0018 

78 

2400 

4 

0.0018 

89 

120 

5 

0.0018 


Figure 6 shows a comparison between GFSSP and Hardy Cross predicted flowrates. The 
comparison appears reasonable considering the fact that Hardy Cross method assumes a 
constant friction factor in the branch while GFSSP computes the friction factor for each 
branch during every iteration. Therefore, as the flowrates change the friction factor also 
changes. 
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Figure 6. - A Flow Rate Comparison Between GFSSP and Hardy Cross Method 
Predictions 

Blowdown of a Tank 

Figure 7 shows a tank with an internal volume of 10 ft^ containing nitrogen gas at a 
pressure and temperature of 100 psia and 80 °F respectively. The nitrogen is discharged 
into the atmosphere through an orifice with a 0.1 inch diameter until the pressure in the 
tank becomes 50 psia. 
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Figure 7. - Physical Schematic (a) and GFSSP Model (b) for Venting Nitrogen Tank 



This is an initial value problem and the initial conditions are: 
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The analytical solution for p / Pj is given by Moody [23] as: 
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The analytical and GFSSP solutions are compared in Figure 8. Excellent agreement was 
observed between the two solutions. 
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Figure 8 - Comparison of the Predicted Pressure History by GFSSP and the 

Analytical Solution 
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